#------------------------------------------------------------------------------
# Import Libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, 
       zoo, tictoc, fst, fixest, 
       PanelMatch, patchwork,
       rio, magrittr, janitor, did, 
       panelView, ggiplot, 
       tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------


#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Regressions
#------------------------------------------------------------------------------

# %%
if (exists("vcf_data")){
  vcf = vcf_data[year>= 1995]
}
vcf[, t := year - 1995]
vcf[, time := year - first_pesa_exposure]


# %% above median sample - cutoff in 1990
evsamp = vcf[cover_1990 > quantile(vcf$cover_1990, 0.5)][time 
                                                        %between% c(-6, 10)]
es1 = feols(forest_index  ~ i(time, sch,  ref=-1) | cellid + styear,
            cluster = ~ blk, evsamp)
es2 = feols(green_index ~ i(time, sch,  ref=-1) | cellid + styear,
            cluster = ~ blk, evsamp)
es3 = feols(built_index ~ i(time, sch,  ref=-1) | cellid + styear,
            cluster = ~ blk, evsamp)

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Generate Tables
#------------------------------------------------------------------------------

# Mean of Y beforetreatment for never treated and treated
evsamp[, never_treated := max(D) == 0, cellid]
controls_pre1 = evsamp[never_treated == 1 & year < first_pesa_exposure, 
                       mean(forest_index)]
treat_pre1    = evsamp[never_treated == 0 & year < first_pesa_exposure, 
                       mean(forest_index)]

# Change name of coefficients
treatmap =c(
  "time::-6:sch" = "Lag Year 6 $\\times$ Sch",
  "time::-5:sch" = "Lag Year 5 $\\times$ Sch",
  "time::-4:sch" = "Lag Year 4 $\\times$ Sch",
  "time::-3:sch" = "Lag Year 3 $\\times$ Sch",
  "time::-2:sch" = "Lag Year 2 $\\times$ Sch",
  "time::0:sch" = "Year 0 $\\times$ Sch",
  "time::1:sch" = "Lead Year 1 $\\times$ Sch",
  "time::2:sch" = "Lead Year 2 $\\times$ Sch",
  "time::3:sch" = "Lead Year 3 $\\times$ Sch",
  "time::4:sch" = "Lead Year 4 $\\times$ Sch",
  "time::5:sch" = "Lead Year 5 $\\times$ Sch",
  "time::6:sch" = "Lead Year 6 $\\times$ Sch",
  "time::7:sch" = "Lead Year 7 $\\times$ Sch",
  "time::8:sch" = "Lead Year 8 $\\times$ Sch",
  "time::9:sch" = "Lead Year 9 $\\times$ Sch",
  "time::10:sch" = "Lead Year 10 $\\times$ Sch",
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-Forest Vegetation",
  "built_index"  = "Bare Ground Indices",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)


fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
desc = "Dynamic Treatment Effects of PESA Adoption on 
Forest Index, Non-Forest Vegetation,
and Bare Ground Indices"
fn = "figure5_regs"; lab = "tab:figure5_regs";
etable( list(es1, es2, es3),
        style.tex = style.tex(main = "base", 
                              depvar.title = "", 
                              model.title = "", 
                              var.title = "\\midrule", 
                              slopes.title = "\\midrule \\emph{Time Trends}", 
                              yesNo = c("$\\checkmark$", ""), 
                              signif.code = NA, 
                              tablefoot = FALSE, 
                              line.bottom = "\\midrule \\midrule"),
        signif.code = NA,
        fixef_sizes = T, 
        fixef_sizes.simplify = F,
        fitstat = c( "n_new", "r2" ), 
        tex = TRUE,
        dict = treatmap,
       label = lab,
       notes = "\\emph{Notes:} Standard errors are clustered at the block level and reported in parentheses.",
       title = glue::glue("{desc}"),
       extralines=list(
         "\\midrule \\emph{Summary Statistics}" = c( "", "", "" ),
         "Mean Pre-Y (Non-Sch)"= c( rep( controls_pre1, 3) ),
         "Mean Pre-Y (Sch )"   = c( rep( treat_pre1, 3 ) ),
         "Dataset"     = c( rep( "VCF", 3 ) ),
         "Timespan"    = c( rep( "1995-2017", 3 ) )
       ),
       file = "appendix_figureA4_table.tex", 
       replace = TRUE
)

#------------------------------------------------------------------------------

